function make_clim_fromROMS(roms_dir, out_dir, newgrid, newS, new_timespan, new_timestep)
%**************************************************
% make_clim_fromROMS.m  
%
% edited by DAS, 7/15/2009 from make_clim.m
%
% This creates a NetCDF file to be used for ROMS forcing.  It specifically reads in 
% history files created by another ROMS run and makes a ROMS input climatology file
% with information on T, S, u, v, ubar, vbar and zeta. It extrapolates by nearest 
% neighbor to fill out the entire field (i.e. it does the preprocess_clim.m step too)
%
% INPUT: 1) roms_dir = directory where ocean_his files from old ROMS run are
%        2) out_dir = where you want your new climatology file to go
%        3) newgrid = name of .nc gridfile for your new grid, generated by make_grid.m
%        4) newS = structure S saved in info.mat for your new run, generated in make_a_run.m
%        5) new_timespan = optional, can specify [start_time end_time] 2-element vector in 
%           MATLAB datenum format for a new time range to write out your climatology 
%        6) new_timestep = optional, though if you specify a new_timespan, you must also 
%           have this variable. it's the number of steps between reading in each history
%           file, so choose 24 if you want to write a new ocean_clm that has daily fields, etc.
%
% NOTE: This file is best run separately from make_a_run.m
%
%**************************************************

new_cname = 'ocean_clm4.nc';  % name of new clm file to make
clmname = [out_dir, new_cname]; %name of new climatology file

% make the directory out_dir if needed
if ~exist(out_dir,'dir'); 
    mkdir(out_dir); 
end

tic

% Read in the original grid matrices
files = dir([roms_dir, 'ocean_his_*']); %get all history files
startfile = [roms_dir, files(1).name];

% if want same time span as original roms run, don't need new_timespan argument
if(nargin < 5); 
    new_timestep = 1; %keep hourly as the history files probably are
    new_timespan = 1:new_timestep:length(files); %keep files as is
else
    endfile = [roms_dir, files(end).name];
    [Gend,Send,Tend] = Z_get_basic_info(endfile); 
    [G,S,T] = Z_get_basic_info(startfile);
    tspan_old = T.time_datenum:1/24:Tend.time_datenum; %assume saved every 1 hour
    tind1 = dsearchn(tspan_old(:), new_timespan(1));
    tind2 = dsearchn(tspan_old(:), new_timespan(2));
    new_timespan = tind1:new_timestep:tind2;
end

% make new file list based on time span you want
for i=1:length(new_timespan)
    newfiles(i).name = files(new_timespan(i)).name;
    vartime(i) = nc_varget([roms_dir,newfiles(i).name],'ocean_time'); %define time vector
end
startfile = [roms_dir, newfiles(1).name];
[G,S,T] = Z_get_basic_info(startfile); %gets grid (G), S, and T from first file

% Read in the new grid matrices
lon_rho = nc_varget(newgrid,'lon_rho'); lat_rho = nc_varget(newgrid,'lat_rho');
lon_u = nc_varget(newgrid,'lon_u'); lat_u = nc_varget(newgrid,'lat_u');
lon_v = nc_varget(newgrid,'lon_v'); lat_v = nc_varget(newgrid,'lat_v');
h = nc_varget(newgrid,'h'); % depth (m) positive down
[Mrho,Lrho] = size(lon_rho);
[Mu,Lu] = size(lon_u);
[Mv,Lv] = size(lon_v);

% set up list of variable names to create climatology files for
varname_list = { ...
    'salt'; ...
    'temp'; ...
    'u'; ...
    'v'; ...
    'zeta'; ...
    'ubar'; ...
    'vbar'; ...
    };
nvar = length(varname_list);

% create the NetCDF output file
my_mode = bitor ( nc_clobber_mode, nc_64bit_offset_mode );
nc_create_empty( clmname, my_mode );
% global attributes
nc_padheader (clmname, 20000 );
nc_attput(clmname, nc_global, 'type','ROMS Climatology File');
% define some dimensions
nc_add_dimension(clmname, 'xi_rho', Lrho); nc_add_dimension(clmname, 'eta_rho', Mrho);
nc_add_dimension(clmname, 'xi_u', Lu); nc_add_dimension(clmname, 'eta_u', Mu);
nc_add_dimension(clmname, 'xi_v', Lv); nc_add_dimension(clmname, 'eta_v', Mv);
nc_add_dimension(clmname, 's_rho', newS.N);
% and we will define the time dimension separately for each variable
ntimes = length(vartime); %number of time steps
  
for ivv = 1:nvar % start of VARNAME LOOP
    do_uvbar = 0; % switch to write uvbar [=1]
    varname = char(varname_list(ivv));
    switch varname
        case 'salt'
            nclongname = 'salinity climatology';
            ncunits = 'PSU';
            nctimename = 'salt_time';
            do_depth = 1; % flag to indicate that there is depth
            do_addtimevar = 1;
            ncxiname = 'xi_rho'; ncetaname = 'eta_rho';
            AAlon = lon_rho; AAlat = lat_rho;
            AAlon_old = G.lon_rho; AAlat_old = G.lat_rho; mask_old = G.mask_rho;
        case 'temp'
            nclongname = 'potential temperature climatology';
            ncunits = 'Celsius';
            nctimename = 'temp_time';
            do_depth = 1; % flag to indicate that there is depth
            do_addtimevar = 1;
            ncxiname = 'xi_rho'; ncetaname = 'eta_rho';
            AAlon = lon_rho; AAlat = lat_rho;
            AAlon_old = G.lon_rho; AAlat_old = G.lat_rho; mask_old = G.mask_rho;
        case 'u'
            nclongname = 'u-momentum component climatology';
            ncunits = 'meter second-1';
            nctimename = 'v3d_time';
            do_depth = 1; % flag to indicate that there is depth
            do_addtimevar = 1;
            ncxiname = 'xi_u'; ncetaname = 'eta_u';
            AAlon = lon_u; AAlat = lat_u;
            AAlon_old = G.lon_u; AAlat_old = G.lat_u; mask_old = G.mask_u;
        case 'v'
            nclongname = 'v-momentum component climatology';
            ncunits = 'meter second-1';
            nctimename = 'v3d_time';
            do_depth = 1; % flag to indicate that there is depth
            do_addtimevar = 0;
            ncxiname = 'xi_v'; ncetaname = 'eta_v';
            AAlon = lon_v; AAlat = lat_v;
            AAlon_old = G.lon_v; AAlat_old = G.lat_v; mask_old = G.mask_v;
        case 'zeta'
            nclongname = 'sea surface height climatology';
            ncunits = 'meter';
            nctimename = 'zeta_time';
            do_depth = 0; % flag to indicate that there is depth
            do_addtimevar = 1; 
            ncxiname = 'xi_rho'; ncetaname = 'eta_rho';
            AAlon = lon_rho; AAlat = lat_rho;
            AAlon_old = G.lon_rho; AAlat_old = G.lat_rho; mask_old = G.mask_rho;
        case 'ubar'
            do_uvbar = 1;   % turn this switch on
            nclongname = 'vertically averaged u-momentum climatology';
            ncunits = 'meter second-1';
            nctimename = 'v2d_time';
            do_depth = 0; % flag to indicate that there is depth
            do_addtimevar = 1;
            ncxiname = 'xi_u'; ncetaname = 'eta_u';
            AAlon = lon_u; AAlat = lat_u;
            AAlon_old = G.lon_u; AAlat_old = G.lat_u; mask_old = G.mask_u;
        case 'vbar'
            do_uvbar = 1;   % turn this switch on
            nclongname = 'vertically averaged v-momentum climatology';
            ncunits = 'meter second-1';
            nctimename = 'v2d_time';
            do_depth = 0; % flag to indicate that there is depth
            do_addtimevar = 0;
            ncxiname = 'xi_v'; ncetaname = 'eta_v';
            AAlon = lon_v; AAlat = lat_v;
            AAlon_old = G.lon_v; AAlat_old = G.lat_v; mask_old = G.mask_v;
    end
    [M,L] = size(AAlon); [M_old,L_old] = size(AAlon_old);
    
    disp(['... Writing ',varname]);

    if ~do_uvbar
        %next find indices within time span of interest (vartime is in seconds)
        switch varname
            case 'u'
                uvvartime = vartime;
                ubar_mat = NaN * ones(ntimes,M,L);
            case 'v';
                uvvartime = vartime;
                vbar_mat = NaN * ones(ntimes,M,L);
        end
    end

    if do_uvbar
        vartime = uvvartime;
    end
    
    if(do_addtimevar)
    % make more of NetCDF output file
    % define rest of variables and attributes
     nc_add_dimension(clmname, nctimename, ntimes); %dimension
     varstruct.Name = nctimename;
     varstruct.Dimension = {nctimename};
        long_name = [nclongname,' time'];
        units = 'seconds';
        varstruct.Attribute = struct('Name', ...
            {'long_name','units'},'Value',{long_name,units});
     nc_addvar(clmname, varstruct);
     nc_varput(clmname, nctimename, vartime'); %because now all time var's are unlimited
    end

    %   then the field
    if do_depth
        varstruct.Name = varname;
        varstruct.Dimension = {nctimename, 's_rho', ncetaname, ncxiname};
        long_name = nclongname;
        units = ncunits;
        varstruct.Attribute = struct('Name', ...
            {'long_name','units'},'Value',{long_name,units});
        nc_addvar(clmname, varstruct);
        indepth = nc_varget(startfile, 's_rho'); %get depth info from original grid
        AN = length(indepth);
    else
        varstruct.Name = varname;
        varstruct.Dimension = {nctimename, ncetaname, ncxiname};
        long_name = nclongname;
        units = ncunits;
        varstruct.Attribute = struct('Name', ...
            {'long_name','units'},'Value',{long_name,units});
        nc_addvar(clmname, varstruct);
    end
 
    if ~do_uvbar
        for tt = 1:length(vartime) % start of TIME LOOP
            infile = [roms_dir,newfiles(tt).name];
            disp(['  tt = ',num2str(tt),' out of ', num2str(ntimes)])
            if tt == 1
                % these are the original lat and lon for this variable
                ALON = AAlon_old; 
                ALAT = AAlat_old; 
                % and we make the new ROMS depths on the relevant grid
                hh = interp2(lon_rho, lat_rho, h, AAlon, AAlat);
                [z_rho,z_w] = Z_s2z(hh,0*hh,newS);
                hh_old = interp2(G.lon_rho, G.lat_rho, G.h, AAlon, AAlat);
                [z_rho_old, z_w_old] = Z_s2z(hh_old,0*hh_old,S);
            end

            % reinitialize the storage arrays
            if do_depth
                % AA is the empty matrix on the new ROMS grid in which to put
                % the field from this time step, but it has the vertical grid
                % of the original ROMS file
                AA = NaN * ones(AN,M,L);
                % AAA is the empty matrix on the ROMS grid in which to put
                % the field from this time step, and it has the vertical grid
                % of the ROMS model.
                AAA = NaN * ones(newS.N,M,L);
            else
                AA = NaN * ones(M,L);
                AAA = AA;
            end

            % load the variable field at this time step
            clear A AAnew this_A
            % INTERPOLATE HORIZONTALLY TO ROMS GRIDS (fills whole domain,
            % if there is ANY data on that level)
            if do_depth
                A = nc_varget(infile, varname); %read in variable from old roms file
                for zz = 1:length(indepth) % start of Z-LOOP
                    this_A = squeeze(A(zz,:,:));
                    for jj = 1:M_old %now go through each row and get rid of NaN's by nearest neighbor
                            this_mask = logical(mask_old(jj,:));
                            if sum(this_mask)>1
                                AAnew(jj,:) = interp1(AAlon_old(jj,this_mask),this_A(jj,this_mask), ...
                                                        AAlon_old(jj,:), 'nearest','extrap');
                            elseif sum(this_mask)==1
                                AAnew(jj,:) = this_A(jj,this_mask);
                            else
                                AAnew(jj,:) = NaN;
                            end
                    end
                    this_AA = interp2(ALON, ALAT, AAnew, AAlon,AAlat);
                    % extensions to the East
                    switch varname
                        case 'salt'
                            fillv = 0;
                            this_AA(isnan(this_AA)) = fillv;
                        case 'temp'
                           fillv = 20;
                            this_AA(isnan(this_AA)) = fillv;
                        case {'u','v'}
                            fillv = 0;
                            this_AA(isnan(this_AA)) = fillv;
                    end
                    AA(zz,:,:) = this_AA;
                end
            else
                A = nc_varget(infile, varname);
                for jj = 1:M_old %now go through each row and get rid of NaN's by nearest neighbor
                            this_mask = logical(mask_old(jj,:));
                            if sum(this_mask)>1
                                AAnew(jj,:) = interp1(AAlon_old(jj,this_mask),A(jj,this_mask), ...
                                                        AAlon_old(jj,:), 'nearest','extrap');
                            elseif sum(this_mask)==1
                                AAnew(jj,:) = A(jj,this_mask);
                            else
                                AAnew(jj,:) = NaN;
                            end
                end
                this_AA = interp2(ALON, ALAT, AAnew, AAlon, AAlat);
                switch varname
                    case 'zeta'
                        fillv = 0;
                        this_AA(isnan(this_AA)) = fillv;
                end
                AA = this_AA;
            end

            % INTERPOLATE VERTICALLY TO new ROMS S-SURFACES
            if do_depth
                for jj = 1:M % start of XY-LOOP
                    for ii = 1:L
                        this_AA = squeeze(AA(:,jj,ii));
                        this_z_roms_old = squeeze(z_rho_old(:,jj,ii));
                        this_z_roms = squeeze(z_rho(:,jj,ii));
                        % pad the top up to zero depth
                        this_z_roms_old = [this_z_roms_old; 0];
                        this_AA = [this_AA; this_AA(end)];
                        ind = find(this_AA==fillv);
                          if(~isempty(ind) && (min(ind)-1)>0) %fill in bottom values
                             this_AA(ind) = this_AA(min(ind)-1);
                          end
                        %this extraps to bottom for ncom data that is 0 sometimes
                        this_AAA = interp1(this_z_roms_old, this_AA, this_z_roms,'nearest','extrap');
                        AAA(:,jj,ii) = this_AAA;
                    end
                end
            else
                AAA = AA;
            end

            % figure out the ubar or vbar field at this time step, and write it
            % to a matrix
            if tt == 1 && do_depth
                DZ = abs(diff(z_w,1,1));
            end
            switch varname
                case 'u'
                    ubar_mat(tt,:,:) = squeeze(sum(AAA.*DZ,1)) ./ hh;
                    if tt == length(vartime)
                        tempufile = [out_dir 'temporary_ubar_storage.mat'];
                        save(tempufile,'ubar_mat')
                    end
                case 'v'
                    vbar_mat(tt,:,:) = squeeze(sum(AAA.*DZ,1)) ./ hh;
                    if tt == length(vartime)
                        tempvfile = [out_dir 'temporary_vbar_storage.mat'];
                        save(tempvfile,'vbar_mat')
                    end
            end
       
            % write the field at one time level
            if do_depth
                [nz,nx,ny]=size(AAA);
                nc_varput(clmname, varname, AAA, [tt-1 0 0 0], [1 nz nx ny]); 
                clear nz nx ny AAA %new_data
            else
                [nx,ny]=size(AAA);
                nc_varput(clmname, varname, AAA, [tt-1 0 0], [1 nx ny]);
                clear nx ny AAA %new_data
            end

        end % end of TIME LOOP
    end

    % get the field for ubar or vbar, and get out of the time loop
    if do_uvbar
        switch varname
            case 'ubar'
                load(tempufile)
                [nt,nx,ny]=size(ubar_mat);
                nc_varput(clmname, varname, ubar_mat, [0 0 0], [length(vartime) nx ny]);
            case 'vbar'
                load(tempvfile)
                [nt,nx,ny]=size(vbar_mat);
                nc_varput(clmname, varname, vbar_mat, [0 0 0], [length(vartime) nx ny]);
        end
    end

end % end of VARNAME LOOP


disp('DONE')

toc

